home *** CD-ROM | disk | FTP | other *** search
/ Aminet 1 (Walnut Creek) / Aminet - June 1993 [Walnut Creek].iso / aminet / misc / sci / ephem_src_4_28.lha / circum.c < prev    next >
C/C++ Source or Header  |  1992-04-17  |  11KB  |  374 lines
  1. /* fill in a Sky struct with all we know about each object.
  2.  *(the user defined objects are in obj.c)
  3.  */
  4.  
  5. #include <stdio.h>
  6. #include <math.h>
  7. #include "astro.h"
  8. #include "circum.h"
  9. #include "screen.h"    /* just for SUN and MOON */
  10.  
  11. /* find body p's circumstances now.
  12.  * to save some time the caller may specify a desired accuracy, in arc seconds.
  13.  * if, based on its mean motion, it would not have moved this much since the
  14.  * last time we were called we only recompute altitude and azimuth and avoid
  15.  * recomputing the planet's heliocentric position. use 0.0 for best possible.
  16.  * we always recompute the user-defined objects' position regardless.
  17.  * return 0 if only alt/az changes, else 1 if all other stuff updated too.
  18.  * N.B: values are for opposition, ie, at fastest retrograde.
  19.  */
  20. body_cir (p, as, np, sp)
  21. int p;
  22. double as;
  23. Now *np;
  24. Sky *sp;
  25. {
  26.     typedef struct {
  27.         double l_dpas;    /* mean days per arc second */
  28.         Now l_now;        /* when l_sky was found */
  29.         double l_ra, l_dec;    /* the eod, ie, unprecessed, ra/dec values */
  30.         Sky l_sky;
  31.     } Last;
  32.     /* must be in same order as the astro.h object #define's */
  33.     static Last last[8] = {
  34.         {.000068, {NOMJD}},    /* mercury */
  35.         {.00017, {NOMJD}},    /* venus */
  36.         {.00015, {NOMJD}},    /* mars */
  37.         {.0012, {NOMJD}},    /* jupiter */
  38.         {.0024, {NOMJD}},    /* saturn */
  39.         {.0051, {NOMJD}},    /* uranus */
  40.         {.0081, {NOMJD}},    /* neptune */
  41.         {.011, {NOMJD}}    /* pluto */
  42.     };
  43.     Last objxlast, objylast;
  44.     double lst, alt, az;
  45.     double ehp, ha, dec;    /* ehp: angular dia of earth from body */
  46.     Last *lp;
  47.     int new;
  48.  
  49.     switch (p) {
  50.     case SUN: return (sun_cir (as, np, sp));
  51.     case MOON: return (moon_cir (as, np, sp));
  52.     case OBJX: lp = &objxlast; break;
  53.     case OBJY: lp = &objylast; break;
  54.     default: lp = last + p; break;
  55.     }
  56.  
  57.     /* if less than l_every days from last time for this planet
  58.      * just redo alt/az.
  59.      * ALWAYS redo objects x and y.
  60.      */
  61.     if (p != OBJX && p != OBJY && same_cir (np, &lp->l_now)
  62.               && about_now (np, &lp->l_now, as*lp->l_dpas)) {
  63.         *sp = lp->l_sky;
  64.         new = 0;
  65.     } else {
  66.         double lpd0, psi0;    /* heliocentric ecliptic long and lat */
  67.         double rp0;        /* dist from sun */
  68.         double rho0;    /* dist from earth */
  69.         double lam, bet;    /* geocentric ecliptic long and lat */
  70.         double dia, mag;    /* angular diameter at 1 AU and magnitude */
  71.         double lsn, rsn;    /* true geoc lng of sun, dist from sn to earth*/
  72.         double el;    /* elongation */
  73.         double f;   /* phase from earth */
  74.  
  75.         lp->l_now = *np;
  76.         sunpos (mjd, &lsn, &rsn);
  77.         if (p == OBJX || p == OBJY)
  78.         obj_cir(mjd, p, &lpd0, &psi0, &rp0, &rho0, &lam, &bet,
  79.                         &sp->s_size, &sp->s_mag);
  80.         else {
  81.         double deps, dpsi;
  82.         double a;
  83.         plans(mjd, p, &lpd0, &psi0, &rp0, &rho0, &lam, &bet, &dia,&mag);
  84.         nutation (mjd, &deps, &dpsi);    /* correct for nutation */
  85.         lam += dpsi;
  86.         a = lsn-lam;            /* and 20.4" aberation */
  87.         lam -= degrad(20.4/3600)*cos(a)/cos(bet);
  88.         bet -= degrad(20.4/3600)*sin(a)*sin(bet);
  89.         }
  90.  
  91.         ecl_eq (mjd, bet, lam, &lp->l_ra, &lp->l_dec);
  92.  
  93.         sp->s_ra = lp->l_ra;
  94.         sp->s_dec = lp->l_dec;
  95.         if (epoch != EOD)
  96.         precess (mjd, epoch, &sp->s_ra, &sp->s_dec);
  97.         sp->s_edist = rho0;
  98.         sp->s_sdist = rp0;
  99.         elongation (lam, bet, lsn, &el);
  100.         el = raddeg(el);
  101.         sp->s_elong = el;
  102.         f = (rp0 > 0.0)
  103.         ? 0.25 * (((rp0+rho0)*(rp0+rho0) - rsn*rsn)/(rp0*rho0)) : 0.0;
  104.         sp->s_phase = f*100.0; /* percent */
  105.         if (p != OBJX && p != OBJY) {
  106.         sp->s_size = dia/rho0;
  107.         sp->s_mag = mag + 5.0*log(rp0*rho0/sqrt(f))/log(10.0);
  108.         }
  109.         sp->s_hlong = lpd0;
  110.         sp->s_hlat = psi0;
  111.         new = 1;
  112.     }
  113.  
  114.     /* alt, az; correct for parallax and refraction; use eod ra/dec */
  115.     now_lst (np, &lst);
  116.     ha = hrrad(lst) - lp->l_ra;
  117.     if (sp->s_edist > 0.0) {
  118.         ehp = (2.0*6378.0/146.0e6) / sp->s_edist;
  119.         ta_par (ha, lp->l_dec, lat, height, ehp, &ha, &dec);
  120.     } else
  121.         dec = lp->l_dec;
  122.     hadec_aa (lat, ha, dec, &alt, &az);
  123.     refract (pressure, temp, alt, &alt);
  124.     sp->s_alt = alt;
  125.     sp->s_az = az;
  126.     lp->l_sky = *sp;
  127.     return (new);
  128. }
  129.  
  130. /* find local times when sun is 18 degrees below horizon.
  131.  * return 0 if just returned same stuff as previous call, else 1 if new.
  132.  */
  133. twilight_cir (np, dawn, dusk, status)
  134. Now *np;
  135. double *dawn, *dusk;
  136. int *status;
  137. {
  138.     static Now last_now = {NOMJD};
  139.     static double last_dawn, last_dusk;
  140.     static int last_status;
  141.     int new;
  142.  
  143.     if (same_cir (np, &last_now) && same_lday (np, &last_now)) {
  144.         *dawn = last_dawn;
  145.         *dusk = last_dusk;
  146.         *status = last_status;
  147.         new = 0;
  148.     } else {
  149.         double x;
  150.         (void) riset_cir (SUN,np,0,TWILIGHT,dawn,dusk,&x,&x,&x,&x,status);
  151.         last_dawn = *dawn;
  152.         last_dusk = *dusk;
  153.         last_status = *status;
  154.         last_now = *np;
  155.         new = 1;
  156.     }
  157.     return (new);
  158. }
  159.  
  160. /* find sun's circumstances now.
  161.  * as is the desired accuracy, in arc seconds; use 0.0 for best possible.
  162.  * return 0 if only alt/az changes, else 1 if all other stuff updated too.
  163.  */
  164. sun_cir (as, np, sp)
  165. double as;
  166. Now *np;
  167. Sky *sp;
  168. {
  169.     static Sky last_sky;
  170.     static Now last_now = {NOMJD};
  171.     static double last_ra, last_dec;    /* unprecessed ra/dec */
  172.     double lst, alt, az;
  173.     double ehp, ha, dec;    /* ehp: angular dia of earth from body */
  174.     int new;
  175.  
  176.     if (same_cir (np, &last_now) && about_now (np, &last_now, as*.00028)) {
  177.         *sp = last_sky;
  178.         new = 0;
  179.     } else {
  180.         double lsn, rsn;
  181.         double deps, dpsi;
  182.  
  183.         last_now = *np;
  184.         sunpos (mjd, &lsn, &rsn);        /* sun's true ecliptic long
  185.                          * and dist
  186.                          */
  187.         nutation (mjd, &deps, &dpsi);    /* correct for nutation */
  188.         lsn += dpsi;
  189.         lsn -= degrad(20.4/3600);        /* and light travel time */
  190.  
  191.         sp->s_edist = rsn;
  192.         sp->s_sdist = 0.0;
  193.         sp->s_elong = 0.0;
  194.         sp->s_size = raddeg(4.65242e-3/rsn)*3600*2;
  195.         sp->s_mag = -26.8;
  196.         sp->s_hlong = lsn-PI;    /* geo- to helio- centric */
  197.         range (&sp->s_hlong, 2*PI);
  198.         sp->s_hlat = 0.0;
  199.  
  200.         ecl_eq (mjd, 0.0, lsn, &last_ra, &last_dec);
  201.         sp->s_ra = last_ra;
  202.         sp->s_dec = last_dec;
  203.         if (epoch != EOD)
  204.         precess (mjd, epoch, &sp->s_ra, &sp->s_dec);
  205.         new = 1;
  206.     }
  207.  
  208.     now_lst (np, &lst);
  209.     ha = hrrad(lst) - last_ra;
  210.     ehp = (2.0 * 6378.0 / 146.0e6) / sp->s_edist;
  211.     ta_par (ha, last_dec, lat, height, ehp, &ha, &dec);
  212.     hadec_aa (lat, ha, dec, &alt, &az);
  213.     refract (pressure, temp, alt, &alt);
  214.     sp->s_alt = alt;
  215.     sp->s_az = az;
  216.     last_sky = *sp;
  217.     return (new);
  218. }
  219.  
  220. /* find moon's circumstances now.
  221.  * as is the desired accuracy, in arc seconds; use 0.0 for best possible.
  222.  * return 0 if only alt/az changes, else 1 if all other stuff updated too.
  223.  */
  224. moon_cir (as, np, sp)
  225. double as;
  226. Now *np;
  227. Sky *sp;
  228. {
  229.     static Sky last_sky;
  230.     static Now last_now = {NOMJD};
  231.     static double ehp;
  232.     static double last_ra, last_dec;    /* unprecessed */
  233.     double lst, alt, az;
  234.     double ha, dec;
  235.     int new;
  236.  
  237.     if (same_cir (np, &last_now) && about_now (np, &last_now, as*.000021)) {
  238.         *sp = last_sky;
  239.         new = 0;
  240.     } else {
  241.         double lam, bet;
  242.         double deps, dpsi;
  243.         double lsn, rsn;    /* sun long in rads, earth-sun dist in au */
  244.         double edistau;    /* earth-moon dist, in au */
  245.         double el;        /* elongation, rads east */
  246.  
  247.         last_now = *np;
  248.         moon (mjd, &lam, &bet, &ehp);    /* moon's true ecliptic loc */
  249.         nutation (mjd, &deps, &dpsi);    /* correct for nutation */
  250.         lam += dpsi;
  251.         range (&lam, 2*PI);
  252.  
  253.         sp->s_edist = 6378.14/sin(ehp);    /* earth-moon dist, want km */
  254.         sp->s_size = 3600*31.22512*sin(ehp);/* moon angular dia, seconds */
  255.         sp->s_hlong = lam;            /* save geo in helio fields */
  256.         sp->s_hlat = bet;
  257.  
  258.         ecl_eq (mjd, bet, lam, &last_ra, &last_dec);
  259.         sp->s_ra = last_ra;
  260.         sp->s_dec = last_dec;
  261.         if (epoch != EOD)
  262.         precess (mjd, epoch, &sp->s_ra, &sp->s_dec);
  263.  
  264.         sunpos (mjd, &lsn, &rsn);
  265.         range (&lsn, 2*PI);
  266.         elongation (lam, bet, lsn, &el);
  267.  
  268.         /* solve triangle of earth, sun, and elongation for moon-sun dist */
  269.         edistau = sp->s_edist/1.495979e8; /* km -> au */
  270.         sp->s_sdist =
  271.         sqrt (edistau*edistau + rsn*rsn - 2.0*edistau*rsn*cos(el));
  272.  
  273.         /* TODO: improve mag; this is based on a flat moon model. */
  274.         sp->s_mag = -12.7 + 2.5*(log10(PI) - log10(PI/2*(1+1.e-6-cos(el))));
  275.  
  276.         sp->s_elong = raddeg(el);    /* want degrees */
  277.         sp->s_phase = fabs(el)/PI*100.0;    /* want non-negative % */
  278.         new = 1;
  279.     }
  280.  
  281.     /* show topocentric alt/az by correcting ra/dec for parallax 
  282.      * as well as refraction.
  283.      */
  284.     now_lst (np, &lst);
  285.     ha = hrrad(lst) - last_ra;
  286.     ta_par (ha, last_dec, lat, height, ehp, &ha, &dec);
  287.     hadec_aa (lat, ha, dec, &alt, &az);
  288.     refract (pressure, temp, alt, &alt);
  289.     sp->s_alt = alt;
  290.     sp->s_az = az;
  291.     last_sky = *sp;
  292.     return (new);
  293. }
  294.  
  295. /* given geocentric ecliptic longitude and latitude, lam and bet, of some object
  296.  * and the longitude of the sun, lsn, find the elongation, el. this is the
  297.  * actual angular separation of the object from the sun, not just the difference
  298.  * in the longitude. the sign, however, IS set simply as a test on longitude
  299.  * such that el will be >0 for an evening object <0 for a morning object.
  300.  * to understand the test for el sign, draw a graph with lam going from 0-2*PI
  301.  *   down the vertical axis, lsn going from 0-2*PI across the hor axis. then
  302.  *   define the diagonal regions bounded by the lines lam=lsn+PI, lam=lsn and
  303.  *   lam=lsn-PI. the "morning" regions are any values to the lower left of the
  304.  *   first line and bounded within the second pair of lines.
  305.  * all angles in radians.
  306.  */
  307. elongation (lam, bet, lsn, el)
  308. double lam, bet, lsn;
  309. double *el;
  310. {
  311.     *el = acos(cos(bet)*cos(lam-lsn));
  312.     if (lam>lsn+PI || lam>lsn-PI && lam<lsn) *el = - *el;
  313. }
  314.  
  315. /* return whether the two Nows are for the same observing circumstances. */
  316. same_cir (n1, n2)
  317. register Now *n1, *n2;
  318. {
  319.     return (n1->n_lat == n2->n_lat
  320.         && n1->n_lng == n2->n_lng
  321.         && n1->n_temp == n2->n_temp
  322.         && n1->n_pressure == n2->n_pressure
  323.         && n1->n_height == n2->n_height
  324.         && n1->n_tz == n2->n_tz
  325.         && n1->n_epoch == n2->n_epoch);
  326. }
  327.  
  328. /* return whether the two Nows are for the same LOCAL day */
  329. same_lday (n1, n2)
  330. Now *n1, *n2;
  331. {
  332.     return (mjd_day(n1->n_mjd - n1->n_tz/24.0) ==
  333.                mjd_day(n2->n_mjd - n2->n_tz/24.0)); 
  334. }
  335.  
  336. /* return whether the mjd of the two Nows are within dt */
  337. static
  338. about_now (n1, n2, dt)
  339. Now *n1, *n2;
  340. double dt;
  341. {
  342.     return (fabs (n1->n_mjd - n2->n_mjd) <= dt/2.0);
  343. }
  344.  
  345. now_lst (np, lst)
  346. Now *np;
  347. double *lst;
  348. {
  349.     utc_gst (mjd_day(mjd), mjd_hr(mjd), lst);
  350.     *lst += radhr(lng);
  351.     range (lst, 24.0);
  352. }
  353.  
  354. /* round a time in days, *t, to the nearest second, IN PLACE. */
  355. rnd_second (t)
  356. double *t;
  357. {
  358.     *t = floor(*t*SPD+0.5)/SPD;
  359. }
  360.     
  361. double
  362. mjd_day(jd)
  363. double jd;
  364. {
  365.     return (floor(jd-0.5)+0.5);
  366. }
  367.  
  368. double
  369. mjd_hr(jd)
  370. double jd;
  371. {
  372.     return ((jd-mjd_day(jd))*24.0);
  373. }
  374.